Women in Work Force

Beyond the Gender Wage Gap

Much attention, both academically and journalistically, has been paid to the “Gender Wage Gap” (e.g. http://inequality.stanford.edu/_media/pdf/key_issues/gender_research.pdf). As this is well-tred territory, our team is interested in looking beyond the gender wage gap and analyzing the inter-gender wage gap. That is, within women in the workforce, are there identifiables gaps in earnings? The following report both confirms a gender wage gap across occupations and identifies gaps among several characteristics of women working within each occupation.

Where is the gap?

Since we know there is a significant wage gap between different professions, we used occupation as a category for our data. We seperated the occupations by groups. Within each group, we divided the men and women and compared their weighted mean.

We started by importing the data set from the United States Census Bureau.

#Importing Data Set
# www2.census.gov/programs-surveys/acs/data/pums/2014/1-Year/csv_pus.zip
ss14pusa = read.csv(file.choose(),header=TRUE)  # read csv file
ss14pusb = read.csv(file.choose(),header=TRUE)

For our anaylsis, we used dplyr, plyr, dplyr, rcokeh, ggplot2, plotly, car, survey, and googleVis.

The following shows the columns we imported.

attach(ss14pusa)
#names(ss14pusa)
ss14pusa_edit = data.frame(SEX,OCCP,WKHP,MAR,WAGP,SCHL,ESP, RAC1P, ST, PAOC, PWGTP, AGEP)
detach(ss14pusa)

attach(ss14pusb)
ss14pusb_edit = data.frame(SEX,OCCP,WKHP,MAR,WAGP,SCHL,ESP, RAC1P, ST, PAOC, PWGTP, AGEP)
detach(ss14pusb)

Data = rbind(ss14pusa_edit,ss14pusb_edit)
colnames(Data) <- c("Gender","Occupation", "Work_hours", "Marriage", "Income", "Education" ,"Parental_Occupation", "Race", "State", "Children", "Weight","Age")

#write.csv(Data, file = "Data.csv",row.names=TRUE)

#Delete income=0/NA rows
row_to_keep=which(Data$Income>0)
Data=Data[row_to_keep,]

We renamed and categorized multiple categories for readablilty.

Note: we used eight high paying occupation to categorize jobs

attach(ss14pusa)
#names(ss14pusa)
ss14pusa_edit = data.frame(SEX,OCCP,WKHP,MAR,WAGP,SCHL,ESP, RAC1P, ST, PAOC, PWGTP, AGEP)
detach(ss14pusa)
attach(ss14pusb)
ss14pusb_edit = data.frame(SEX,OCCP,WKHP,MAR,WAGP,SCHL,ESP, RAC1P, ST, PAOC, PWGTP, AGEP)
detach(ss14pusb)
Data = rbind(ss14pusa_edit,ss14pusb_edit)
colnames(Data) <- c("Gender","Occupation", "Work_hours", "Marriage", "Income", "Education" ,"Parental_Occupation", "Race", "State", "Children", "Weight","Age")
#write.csv(Data, file = "Data.csv",row.names=TRUE)
#Delete income=0/NA rows
row_to_keep=which(Data$Income>0)
Data=Data[row_to_keep,]

Based on the eight occupations, we plot the weighted mean of men and women to compare the gap.

NAs introduced by coercionNAs introduced by coercion

From the graph, we see there is a significant difference between men and women over all occupations. To gain a better understanding we graphed them by percentage from the weighted mean of both genders.

#median wage for each occupation by gender (just to acknowledge that there *is* a gender gap)
Data_sex_occp <- ddply(Data, .(Gender, Occupation), summarise, MEAN = weighted.mean(Income, Weight, na.rm = T))
ggplot(Data_sex_occp, aes(x=Occupation, y=MEAN, fill=factor(Gender))) + 
  geom_bar(stat="identity",position="dodge") + 
  scale_fill_brewer(palette="RdYlGn") +
  labs(fill="") + 
  ylab("Mean Salary ($)") + 
  xlab("Occupations") + 
  ggtitle(paste("Salary Comparison between Men & Women")) + 
  theme(axis.text.x = element_text(angle = 30, hjust = 1), 
        panel.background = element_rect(fill = 'white' )) + 
  theme_grey(base_size = 12)

Oberservations:

1. Medical occupations has the largest disperity between men and women

2. Women in legal occupations are the closest to the median

In this report, we will investigate 5 probable causes to determine which groups of women will most likely succeed or fail. These categories are: Marriage Status, Motherhood, Location, Race, and Education.

Familial Responsibilities

Our next variable of interest is women’s marital statuses. A common explanation of the gender wage gap is that women make the decision with their families to shift their focus to household responsibilities. We’ll be looking at the effect of marriage and motherhood seperately first, then looking how marital statuses of mothers affect their earnings.

Marital Status

Below we analyse the effect of marital status on mean income for women across all occupations, as well as within occupations.

Observations:

1. Overall,widowed and seperated women make significantly less than married women when occupations are pooled.

2. Marital status has no effect on earning for women in the scientific field

3. Interestingly, for most occupations where marital status is significant, women in their profession who have never been married earn more than married women, whereas widowed and seperated women earn significantly less. Perhaps widowed or seperated women are more likely to have children than never married women, and thus must shoulder more responsibilities outside of work.

Motherhood

After investigating how marriage status affect women’s income, we looked into the presence and age of children owned by the women to determine how it might correlate to women’s income. We expect women with Children under 6 may have less working hours and thus lower income (reasoning: when a child is young the women of the house is more likely to stay home with their children. As we see the children grow, we expect a trend for women to earn the higher.) Furthermore, we may also want to figure out how occupation and age for women with children affect their income.

At first, we want to have a glance at the income level of the different groups of females. Basically females are divided into four categories regarding the presence and age of Children: Females with Children under 6, Females with Children from 6 to 17, Females with Children under 6 and above 6, Females with No Children. Thus we simply creat a bar chart to illustrate the income difference.

#marital status, subset at specific age to control for experience
Data_women <- filter(Data, Gender == 2)
Data_wed_occp <- ddply(Data_women, .(Marriage, Occupation, Age), summarise, MEAN = weighted.mean(Income, Weight, na.rm = T))
plot(Data_wed_occp$Age, Data_wed_occp$MEAN)

summary(lm(MEAN ~ Marriage + Age + I(Age^2) , data = Data_wed_occp))

Call:
lm(formula = MEAN ~ Marriage + Age + I(Age^2), data = Data_wed_occp)

Residuals:
   Min     1Q Median     3Q    Max 
-56107 -11573  -1470   8035 231989 

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)    
(Intercept)           -39977.664   3357.450 -11.907  < 2e-16 ***
MarriageWidowed        -5058.692   1430.698  -3.536 0.000415 ***
MarriageDivorced       -3058.190   1342.701  -2.278 0.022842 *  
MarriageSeperated      -8106.266   1479.165  -5.480 4.72e-08 ***
MarriageNever Married    557.144   1337.027   0.417 0.676935    
Age                     4089.016    134.708  30.355  < 2e-16 ***
I(Age^2)                 -38.882      1.275 -30.488  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 21170 on 2232 degrees of freedom
Multiple R-squared:  0.2979,    Adjusted R-squared:  0.296 
F-statistic: 157.8 on 6 and 2232 DF,  p-value: < 2.2e-16

#Finance
Women_FIN<-filter(Data_women, Occupation =="FIN")
Data_wed_FIN <- ddply(Women_FIN, .(Marriage, Age), summarise, MEAN = weighted.mean(Income, Weight, na.rm = T))
FIN <- plot_ly(Data_wed_FIN, x = Age, y = MEAN, mode = "markers", opacity = .4, color = Marriage, title = "Women in Finance") 
FIN <-add_trace(FIN,y = fitted(loess(MEAN ~ Age + as.numeric(Marriage))) , x = Age, color = Marriage) 
layout(FIN, title = "Women in Finance")

summary(lm(MEAN ~ Marriage + Age + I(Age^2), data = Data_wed_FIN))

Call:
lm(formula = MEAN ~ Marriage + Age + I(Age^2), data = Data_wed_FIN)

Residuals:
   Min     1Q Median     3Q    Max 
-35099  -8154   -769   6758  85196 

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)    
(Intercept)           -36346.096   6420.514  -5.661 3.61e-08 ***
MarriageWidowed        -7156.318   2631.738  -2.719  0.00694 ** 
MarriageDivorced       -2261.388   2542.871  -0.889  0.37457    
MarriageSeperated     -15267.284   2763.424  -5.525 7.33e-08 ***
MarriageNever Married   2543.383   2558.912   0.994  0.32108    
Age                     3941.532    258.084  15.272  < 2e-16 ***
I(Age^2)                 -39.041      2.435 -16.032  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 14470 on 291 degrees of freedom
Multiple R-squared:  0.5055,    Adjusted R-squared:  0.4953 
F-statistic: 49.57 on 6 and 291 DF,  p-value: < 2.2e-16

#Engineering
Women_ENG<-filter(Data_women, Occupation =="ENG")
Data_wed_ENG <- ddply(Women_ENG, .(Marriage, Age), summarise, MEAN = weighted.mean(Income, Weight, na.rm = T))
#plotly
ENG <- plot_ly(Data_wed_ENG, x = Age, y = MEAN, mode = "markers", opacity = .4, color = Marriage, title = "Women in Finance") 
ENG <-add_trace(ENG,y = fitted(loess(MEAN ~ Age + as.numeric(Marriage))) , x = Age, color = Marriage) 
layout(ENG, title = "Women in Engineering")

lm1 <- lm(MEAN ~ Marriage + Age + I(Age^2), data = Data_wed_ENG)
summary(lm1)

Call:
lm(formula = MEAN ~ Marriage + Age + I(Age^2), data = Data_wed_ENG)

Residuals:
   Min     1Q Median     3Q    Max 
-71787 -13616  -1922   6880 221919 

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)    
(Intercept)           -67401.667  18767.237  -3.591 0.000411 ***
MarriageWidowed        -1410.202   7393.487  -0.191 0.848921    
MarriageDivorced       -4852.957   6009.505  -0.808 0.420283    
MarriageSeperated     -16473.248   7418.879  -2.220 0.027478 *  
MarriageNever Married  -1392.858   5878.545  -0.237 0.812939    
Age                     5917.840    827.959   7.148 1.50e-11 ***
I(Age^2)                 -59.961      8.655  -6.928 5.36e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 30550 on 206 degrees of freedom
Multiple R-squared:  0.2152,    Adjusted R-squared:  0.1924 
F-statistic: 9.416 on 6 and 206 DF,  p-value: 3.912e-09

#Medicine
Women_MED<-filter(Data_women, Occupation =="MED")
Data_wed_MED <- ddply(Women_MED, .(Marriage, Age), summarise, MEAN = weighted.mean(Income, Weight, na.rm = T))
#plotly
MED <- plot_ly(Data_wed_MED, x = Age, y = MEAN, mode = "markers", opacity = .4, color = Marriage, title = "Women in Finance") 
MED <-add_trace(MED,y = fitted(loess(MEAN ~ Age + as.numeric(Marriage))) , x = Age, color = Marriage) 
layout(MED, title = "Women in Finance")

summary(lm(MEAN ~ Marriage + Age + I(Age^2), data = Data_wed_MED))

Call:
lm(formula = MEAN ~ Marriage + Age + I(Age^2), data = Data_wed_MED)

Residuals:
   Min     1Q Median     3Q    Max 
-35762  -7890    447   6423 117974 

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)    
(Intercept)           -19061.330   6542.696  -2.913 0.003825 ** 
MarriageWidowed       -10048.240   2875.561  -3.494 0.000542 ***
MarriageDivorced       -1420.377   2803.367  -0.507 0.612733    
MarriageSeperated      -9171.092   3002.385  -3.055 0.002441 ** 
MarriageNever Married   -419.022   2828.361  -0.148 0.882317    
Age                     3019.839    255.543  11.817  < 2e-16 ***
I(Age^2)                 -28.188      2.342 -12.037  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 16640 on 323 degrees of freedom
Multiple R-squared:  0.3316,    Adjusted R-squared:  0.3192 
F-statistic:  26.7 on 6 and 323 DF,  p-value: < 2.2e-16

#Science
Women_SCI<-filter(Data_women, Occupation =="SCI")
Data_wed_SCI <- ddply(Women_SCI, .(Marriage, Age), summarise, MEAN = weighted.mean(Income, Weight, na.rm = T))
#plotly
SCI <- plot_ly(Data_wed_SCI, x = Age, y = MEAN, mode = "markers", opacity = .4, color = Marriage, title = "Women in Science") 
SCI <-add_trace(SCI,y = fitted(loess(MEAN ~ Age + as.numeric(Marriage))) , x = Age, color = Marriage) 
layout(SCI, title = "Women in Science")

summary(lm(MEAN ~ Marriage + Age + I(Age^2), data = Data_wed_SCI))

Call:
lm(formula = MEAN ~ Marriage + Age + I(Age^2), data = Data_wed_SCI)

Residuals:
   Min     1Q Median     3Q    Max 
-68959 -10504   -932   8858  95689 

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)    
(Intercept)           -77292.909  10726.610  -7.206 7.04e-12 ***
MarriageWidowed        -6390.342   4489.513  -1.423    0.156    
MarriageDivorced       -2419.314   3934.713  -0.615    0.539    
MarriageSeperated       2882.259   4602.021   0.626    0.532    
MarriageNever Married   2465.396   3899.844   0.632    0.528    
Age                     5640.076    444.813  12.680  < 2e-16 ***
I(Age^2)                 -53.246      4.334 -12.286  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 21450 on 246 degrees of freedom
Multiple R-squared:  0.4066,    Adjusted R-squared:  0.3922 
F-statistic:  28.1 on 6 and 246 DF,  p-value: < 2.2e-16

#Managers
Women_MGR<-filter(Data_women, Occupation =="MGR")
Data_wed_MGR <- ddply(Women_MGR, .(Marriage, Age), summarise, MEAN = weighted.mean(Income, Weight, na.rm = T))
MGR <- plot_ly(Data_wed_MGR, x = Age, y = MEAN, mode = "markers", opacity = .4, color = Marriage, title = "Women in Management") 
MGR <-add_trace(MGR,y = fitted(loess(MEAN ~ Age + as.numeric(Marriage))) , x = Age, color = Marriage) 
layout(MGR, title = "Women in Management")

summary(lm(MEAN ~ Marriage + Age + I(Age^2), data = Data_wed_MGR))

Call:
lm(formula = MEAN ~ Marriage + Age + I(Age^2), data = Data_wed_MGR)

Residuals:
   Min     1Q Median     3Q    Max 
-40481  -9893   -582   5931 212526 

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)    
(Intercept)           -44162.963   8226.978  -5.368 1.51e-07 ***
MarriageWidowed        -2612.919   3617.947  -0.722   0.4707    
MarriageDivorced       -4895.286   3620.382  -1.352   0.1773    
MarriageSeperated      -8413.422   3807.248  -2.210   0.0278 *  
MarriageNever Married  -1116.802   3591.031  -0.311   0.7560    
Age                     4377.509    326.101  13.424  < 2e-16 ***
I(Age^2)                 -40.301      3.009 -13.391  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 21380 on 327 degrees of freedom
Multiple R-squared:  0.3585,    Adjusted R-squared:  0.3467 
F-statistic: 30.46 on 6 and 327 DF,  p-value: < 2.2e-16

#Lawyers
Women_LGL<-filter(Data_women, Occupation =="LGL")
Data_wed_LGL <- ddply(Women_LGL, .(Marriage, Age), summarise, MEAN = weighted.mean(Income, Weight, na.rm = T))
LGL <- plot_ly(Data_wed_LGL, x = Age, y = MEAN, mode = "markers", opacity = .4, color = Marriage) 
LGL <-add_trace(LGL,y = fitted(loess(MEAN ~ Age + as.numeric(Marriage))) , x = Age, color = Marriage) 
layout(LGL, title = "Women in Law")

summary(lm(MEAN ~ Marriage + Age + I(Age^2), data = Data_wed_LGL))

Call:
lm(formula = MEAN ~ Marriage + Age + I(Age^2), data = Data_wed_LGL)

Residuals:
   Min     1Q Median     3Q    Max 
-38365  -8400  -2126   6518 159634 

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)    
(Intercept)           -39807.860   9552.619  -4.167 4.26e-05 ***
MarriageWidowed       -11672.166   3875.330  -3.012  0.00286 ** 
MarriageDivorced       -4590.456   3398.000  -1.351  0.17795    
MarriageSeperated      -3691.867   3913.706  -0.943  0.34644    
MarriageNever Married     50.710   3386.857   0.015  0.98807    
Age                     3548.182    390.813   9.079  < 2e-16 ***
I(Age^2)                 -32.067      3.797  -8.446 2.56e-15 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 18510 on 248 degrees of freedom
Multiple R-squared:  0.274, Adjusted R-squared:  0.2565 
F-statistic:  15.6 on 6 and 248 DF,  p-value: 3.508e-15

#Business
Women_BUS<-filter(Data_women, Occupation =="BUS")
Data_wed_BUS <- ddply(Women_BUS, .(Marriage, Age), summarise, MEAN = weighted.mean(Income, Weight, na.rm = T))
BUS <- plot_ly(Data_wed_BUS, x = Age, y = MEAN, mode = "markers", opacity = .4, color = Marriage) 
BUS <-add_trace(BUS,y = fitted(loess(MEAN ~ Age + as.numeric(Marriage))) , x = Age, color = Marriage) 
layout(BUS, title = "Women in Business")

summary(lm(MEAN ~ Marriage + Age + I(Age^2), data = Data_wed_BUS))

Call:
lm(formula = MEAN ~ Marriage + Age + I(Age^2), data = Data_wed_BUS)

Residuals:
   Min     1Q Median     3Q    Max 
-40066 -10052  -1659   6642 173126 

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)    
(Intercept)           -33142.936   8924.930  -3.714 0.000246 ***
MarriageWidowed         2308.734   3704.024   0.623 0.533580    
MarriageDivorced       -1095.901   3625.929  -0.302 0.762688    
MarriageSeperated     -10299.949   3862.286  -2.667 0.008092 ** 
MarriageNever Married   3826.334   3582.032   1.068 0.286327    
Age                     3701.176    355.322  10.416  < 2e-16 ***
I(Age^2)                 -36.125      3.345 -10.800  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 20210 on 287 degrees of freedom
Multiple R-squared:  0.303, Adjusted R-squared:  0.2884 
F-statistic: 20.79 on 6 and 287 DF,  p-value: < 2.2e-16

#Computer
Women_CMM<-filter(Data_women, Occupation =="CMM")
Data_wed_CMM <- ddply(Women_CMM, .(Marriage, Age), summarise, MEAN = weighted.mean(Income, Weight, na.rm = T))
CMM <- plot_ly(Data_wed_CMM, x = Age, y = MEAN, mode = "markers", opacity = .4, color = Marriage) 
CMM <-add_trace(CMM,y = fitted(loess(MEAN ~ Age + as.numeric(Marriage))) , x = Age, color = Marriage) 
layout(CMM, title = "Women in Law")

summary(lm(MEAN ~ Marriage + Age + I(Age^2), data = Data_wed_CMM))

Call:
lm(formula = MEAN ~ Marriage + Age + I(Age^2), data = Data_wed_CMM)

Residuals:
   Min     1Q Median     3Q    Max 
-49094  -9844     78   6958 121180 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)           -47534.01    8781.27  -5.413 1.43e-07 ***
MarriageWidowed        -2742.11    3575.69  -0.767    0.444    
MarriageDivorced       -2754.09    3335.19  -0.826    0.410    
MarriageSeperated      -4313.13    3625.43  -1.190    0.235    
MarriageNever Married  -1096.28    3334.11  -0.329    0.743    
Age                     4792.97     364.29  13.157  < 2e-16 ***
I(Age^2)                 -46.82       3.59 -13.041  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 18000 on 255 degrees of freedom
Multiple R-squared:  0.4061,    Adjusted R-squared:  0.3921 
F-statistic: 29.06 on 6 and 255 DF,  p-value: < 2.2e-16

Observations:

1.Women with children 6 to 17 has obviously higher Average Income as they have more working experience and do not need to spend too much time to take care of young children.

2.Women with Children of both ages have slightly lower income, which is also reasonable as they may not have too much time in work.

The bar chart verifies our guess. Next, in order to better understand why women with children 6 to 17 have higher income and whether it has anything to do with age, we are going to study how the income change for different groups of women through out age. Basically we try to analyze by line chart to see the income trend and compare.

Data_women <- filter(Data, Gender == 2)
Data_women$Agegroup <- cut(as.numeric(Data_women$Age), breaks = c(20,25,30,35,40,45,50,55,60,65))
Data_women <- filter(Data_women, Age > 20 & Age < 66)
### Line Chart with regard to Age
MotherHood=select(Data_women, Income,Children,Agegroup)
#detach(package:plyr)
GroupedMotherHood <-
 MotherHood %>%
 na.omit() %>%
 group_by(Children,Agegroup) %>%
 summarize(
   AvgIncome = mean(Income)
 )
GroupedMotherHood$Children<-factor(GroupedMotherHood$Children, levels = c("1","2","3","4"), c('Children under 6','Children from 6 to 17','Children under 6 and above 6','No Children'))

Error in `$<-.data.frame`(`*tmp*`, "Children", value = integer(0)) : 
  replacement has 0 rows, data has 1

Observations:

1. Before 45-50 years old, there is an abvious increasing trend of income basically for all groups. And after that age, average income begins to decrease or stablize.

2. Surprisingly, for the same age group, women with Children under 6 have higher income. While women with children from 6 and 17 have higher income after 45 years old, which might be the reason this group has higher average income. On the other hand, women with no children doesn’t earn less than other groups at younger age groups, but their income become stable and less than others after 40 years old.

Next, we wanted to see how the working hours related to the income of women in different group. Will women with children from 6 to 17 have more working hours and generate higher income? Thus we want to use bubble chart to integrate working hours, average income, the group count in a bubble chart so that all the information can be gained in the same graph.

#ggplot(GroupedMotherHood, aes(x=Agegroup, y=AvgIncome,fill=Children)) +
#  geom_bar(width=.5, stat="identity")+
#  xlab("Age Group") + ylab("Mean Income") +
#  ggtitle("Income for Motherhood")
ggplot(GroupedMotherHood, aes(x=Agegroup, y=AvgIncome,group=Children,colour=Children)) +
 geom_line()+
 geom_point()+
 #guides(fill=FALSE) +
 xlab("Age Group") + ylab("Mean Income") +
 ggtitle("Income for Motherhood")

Observations:

1. When work_hours smaller than 60Hrs, Avg Income tend to be positively related to work hour, which agrees with our common knowledge. Most women concentrate around 40 working hours.

2. There are generally more females with no Children in this survey as green circles are generally larger in size. Also it can be observed there are much more women with no children working more than 50 hours compared to other groups. On the other hand, there are basically less women with children from both age groups.

3. Females with Children from 6 to 17 often have higher avg income. While for the same working hour level, females with no children are frequently have lower average income.

Next we want to investigate whether the income difference among four grouops of female have any relationship to their occupation. Therefore we use this chart to compare the income for four groups of women in different occupations. Basically, darker color represent higher income and larger size indicate larger group size.

MotherWorking=select(Data_women, Income,Children,Work_hours)
#detach(package:plyr)
GroupedMotherWorking <-
  MotherWorking %>%
  na.omit() %>%
  group_by(Children,Work_hours) %>% 
  summarize(AvgIncome = mean(Income), count=n()) 
GroupedMotherWorking$Children<-factor(GroupedMotherWorking$Children, levels = c("1","2","3","4"),labels= c('Children under 6','Children from 6 to 17','Children under 6 and above 6','No Children'))
# note how size is automatically scaled and added as hover text
plot_ly(GroupedMotherWorking, x = Work_hours, y = AvgIncome,size=sqrt(count), color = Children,mode = "markers",hoverinfo = "text",text = paste("AvgIncome:",AvgIncome,"<br>","Count:",count,"<br>"))

Marital statuses of mothers

Observations: We can see in the summary results of the interaction between marital status and motherhood, single mothers earn significantly less when occupations are pooled and age is controlled.

Education

We expect that educational attainment would affect the women’s income. We may want to figure out which degree is the most important factor to earn more than the median income. Furthermore, we will see which occupation has the largest/smallest “below median group” within each educational attainment.

We divided female into 2 groups (>=$60,000 vs <60,000) using weighted median income for working class. ($60,000) We will see the proportion of these 2 groups for each educational attainment.

MotherMaritalOccupation=select(Data_women, Income,Children,Occupation, Marriage, Age)
MotherMaritalOccupation <- na.omit(MotherMaritalOccupation)
MotherMaritalOccupation$Children <- ifelse(MotherMaritalOccupation$Children == 4, 0, 1)
#detach(package:plyr)
GroupedMotherMaritalOccupation <-
  MotherMaritalOccupation %>%
  na.omit() %>%
  group_by(Children,Occupation, Marriage, Age) %>% 
  summarize(
    AvgIncome = mean(Income),
    count=n()
  ) 
GroupedMotherMarital <-
  MotherMaritalOccupation %>%
  na.omit() %>%
  group_by(Children, Marriage, Age) %>% 
  summarize(
    AvgIncome = mean(Income),
    count=n()
  ) 
ggplot(GroupedMotherMarital,aes(y = AvgIncome, x =Age,colour=as.factor(Children))) +
geom_point() + geom_smooth(method="loess", fill=NA) + facet_wrap(~Marriage, scales= "free")

summary(lm(AvgIncome ~ Age + I(Age^2) + Children*Marriage, data = GroupedMotherMaritalOccupation))

Call:
lm(formula = AvgIncome ~ Age + I(Age^2) + Children * Marriage, 
    data = GroupedMotherMaritalOccupation)

Residuals:
   Min     1Q Median     3Q    Max 
-65888 -10493  -1879   7328 384166 

Coefficients:
                                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)                    -59323.906   5257.629 -11.283  < 2e-16 ***
Age                              4987.877    254.145  19.626  < 2e-16 ***
I(Age^2)                          -47.777      2.929 -16.313  < 2e-16 ***
Children                         2882.840   1628.541   1.770 0.076801 .  
MarriageWidowed                 -6785.321   1847.179  -3.673 0.000244 ***
MarriageDivorced                -2517.617   1618.366  -1.556 0.119902    
MarriageSeperated               -8634.582   1722.751  -5.012 5.71e-07 ***
MarriageNever Married            3243.048   1600.407   2.026 0.042818 *  
Children:MarriageWidowed        -6511.152   2805.557  -2.321 0.020368 *  
Children:MarriageDivorced       -4680.108   2365.606  -1.978 0.047980 *  
Children:MarriageSeperated      -1087.083   2514.920  -0.432 0.665590    
Children:MarriageNever Married -10540.040   2354.884  -4.476 7.91e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 21470 on 2842 degrees of freedom
Multiple R-squared:  0.2619,    Adjusted R-squared:  0.2591 
F-statistic: 91.69 on 11 and 2842 DF,  p-value: < 2.2e-16

Observations:

1. Without college degree, about 80% of them earn less than median income.

2. With master’s or doctor’s degree, about 40% of them earn less than median income.

Next, we use a linear regression.

#Calculate weighted median to split income below/above
Data_w<-svydesign(id = ~1, weights = ~Data$Weight, data = Data)
svyquantile(~Data$Income, Data_w, c(.25,.5,.75))

             0.25   0.5  0.75
Data$Income 35000 60000 92000

##Median income=$60000
#split income below/above
Data$Income_m<-ifelse(Data$Income>=60000,">=$60,000","<$60,000")
Data_women2<-subset(Data,Gender=='2')
attach(Data_women2)

The following objects are masked from Data_women2 (pos = 3):

    Age, Children, Division, Education, Gender, Income, Marriage, Occupation, Parental_Occupation, Race, Region,
    State, Weight, Work_hours

The following objects are masked from Data_women2 (pos = 4):

    Age, Children, Division, Education, Gender, Income, Marriage, Occupation, Parental_Occupation, Race, Region,
    State, Weight, Work_hours

Data_women2<-data.frame(Education,Occupation,Income_m,Weight)
detach(Data_women2)
#get the summary table-education
Data_women_w<-svydesign(id = ~1, weights = ~Data_women2$Weight, data = Data_women2)
t1<-prop.table(svytable(~Income_m+Education, Data_women_w),2)
barplot(t1, main = "Below/Above median income women by Educational attainment",col = c("mistyrose", "lavender"),
xlab = "Educational attainment", names = c("~g8", "g9~nodegree","ass/bach","prof/mas","doc"), 
ylab = "proportion", legend = c("<$60,000", ">=$60,000"), 
args.legend = list(title = "income", x = "topright", cex = 1), ylim = c(0, 1.5))

Observations:

1. factor2 (grade9~college_no_degree) is not significant.

2. Interpretation for factor3: “associate/bachelor’s degree”, versus “~grade8”(baseline), multiplies the odds of being “>=$60000” group by exp(0.94372)

3. We can apply the same interpretation for factor 4,5,6.

4. Within college, master, doctor’s degree, the chance of being in “>=$60000” group increases at higher educational attainment.

5. However, the gap between master’s and doctor’s degree is relatively small.

To futher understand the effect, we are going to break down each education by occupation.

fit_edu<-svyglm(factor(Income_m) ~ factor(Education), design = Data_women_w,family=binomial)

non-integer #successes in a binomial glm!

summary(fit_edu)

Call:
svyglm(formula = factor(Income_m) ~ factor(Education), design = Data_women_w, 
    family = binomial)

Survey design:
svydesign(id = ~1, weights = ~Data_women2$Weight, data = Data_women2)

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)        -1.30179    0.11634 -11.189  < 2e-16 ***
factor(Education)2 -0.06833    0.11711  -0.584     0.56    
factor(Education)3  0.94372    0.11663   8.091 5.94e-16 ***
factor(Education)4  1.78284    0.11708  15.228  < 2e-16 ***
factor(Education)5  2.03576    0.12210  16.673  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1.000005)

Number of Fisher Scoring iterations: 4

Observations:

1. Within “~grade8”, “FIN:finance” has the largest <$60,000 group.

2. “CMM:computer occs” has the smallest <$60,000 group.

Observations:

1. Within “grade9~college_no_degree”, “SCI:science” has the largest <$60,000 group.

2. “CMM:computer occs” has the smallest <$60,000 group.

Observations:

1. Within “associate/bachelor”, “LGL:legal” has the largest <$60,000 group.

2. “CMM: computer occs” has the smallest <$60,000 group.

Observations:

1. Within “Master”, “LGL:legal” has the largest <$60,000 group.

2. “CMM:computer occs” has the smallest <$60,000 group.

Observations:

1. Within “Doctor”, “LGL:legal” has the largest <$60,000 group.

2. “CMM:computer occs” has the smallest <$60,000 group.

Conclusions:

1. Women without college degree is more likely to earn less than median income.

2. (-) At any educational attainment,CMM:computer occs has the smallest <$60,000 group.

3. (+) With the bachelor/master/doctor’s degree LGL:legal has the largest >=$60,000 group.

Regional

You cannot choose where you were born, but you can choose where you live. Different places sends different messages and attracts different people. In this section we are interested to see regional effect on women’s income in the US. The main question here is: in which locations are women better or worse off?

First, we want to look at women’s income boxplots in all 4 regions to better understand women’s income distribution.

The boxplots of women’s income show that despite that all distributions are highly skewed to the right, the Northeast Region has the highest women’s income level of all 4 regions, followed by the West; Midwest has the lowest women’s income level, whereas South is only slightly better off than that. Then we will break down to state level to investigate further on women’s income level.

Income levels can be affected by many things. Cost of living is one of them. For example, In New York, everthing cost a little more than many other parts of the US, so people will only thrive by earning more money. Therefore, in order to remove the effect of cost of living that differs in different states, we calculate the percentage of median of each woman’s income by state and draw boxplot according to that.

From this irisdescent boxplot by state, we can see that New York definitely goes ahead (despite that its cost of living is already high), so does Connecticut with more large extreme values; however California has the highest upper whisker and largest portion of upper middle income; Utah and Wyoming are the worst off.

Then we want to decide locations where women will thrive most financially (relatively to men), so we plot women’s median of percentage income by state in a US map, which shows the percentage of statewise median that half of women has income below. In perfect equal-gender scenarios, the percentage should be 1, which means half of women have income below the exact statewise median.

This map clearly shows that Utah is where women thrive most financially (relatively to men). Half of women in Utah have income below 68.97% of state median, compared to 100% in the perfect equal-gender situation. Surprisingly (at least to us), D.C. tops among all states with a 93%, followed by Nevada 92% and Arizona 90%.

Observation: Women’s income distributions in each state are all highly skewed to the right with lots of large extreme values.

1. The Northeast Region has the highest women’s income level of all 4 regions in the US.

2. After excluding statewise differences e.g. cost of living, women in New York and Connecticut tend to have higher incomes, California has the highest upper whisker and largest portion of upper middle income, and Utah and Wyoming are the worst off.

3. By comparing women’s income medians with state medians, we see that half of Utah’s women have income below 69% of state median, which is the lowest percentage among all states. Woman in D.C., Nevada, and Arizona have income medians closest to state median, with percentage above 90%, whereas other states’ are much lower.

Race

Over the years, there has been an increase awareness on the women wage gap and racial social issues. We wanted to understand the effects of a women’s race could potentially have on a women’s career.

First, we broke down the data of women by race.

# Map women's median of percentage of median by state
PercM <- aggregate(Percentage_of_Median~State,Perc, median)
names(PercM)[1]<-"name"
PM <- merge(PercM, statenames, by="name")
# light grey boundaries
l <- list(color = toRGB("white"), width = 1)
# specify some map projection/options
g <- list(
  scope = 'usa',
  projection = list(type = 'albers usa'),
  showlakes = TRUE,
  lakecolor = toRGB('white')
)
plot_ly(PM, z = Percentage_of_Median, locations = abbr, text=name, type = 'choropleth',
        locationmode = 'USA-states', color = Percentage_of_Median, colors = 'Blues',
        marker = list(line = l), colorbar = list(title = "Median Percentage")) %>%
  layout(title = '2014 Median of Woman Percentage Income by State', geo = g)

Error in seq.default(min(rng), max(rng), length.out = 10) : 
  'from' cannot be NA, NaN or infinite

From the pie chart, we determined that majority of women in the work force are white. To see how this compares to the national average, we compared this chart to the overall race breakdown of the country. Our assumption is that women should make up half of the working force for all races.

Race_freq <- ftable(Data_Women$Race)
Race_freq <- as.data.frame(Race_freq)
plot_ly(Race_freq, labels = Var1, values = Freq, type = "pie",hole = 0.0,name = "Race") %>%
  layout(title = "Percentage of Women by Race")

Observation:

1. From the table, we can gather that Black women are 11% higher than the half of the Black workforce. This means more black women are working than man.

2. All other races have a relatively small difference.

3. Workforce groups with more men than women are White, Asian, and Others.

We also want to look into how women perform in the occupation. We broke up this into three categories:

1. Women who perform under the weighted mean of women

2. Women who perform above the weighted mean of both genders

3. Women who perform above the weighted mean of men

These categories will be classified as below, meet, and exceed respectively.

Race_freq_All <- ftable(Data$Race)
Race_freq_All <- as.data.frame(Race_freq_All)
Race_Compare <- data.frame(Race_freq$Var1,Race_freq$Freq,0.5 *Race_freq_All$Freq,100*(Race_freq$Freq-0.5 *Race_freq_All$Freq)/Race_freq_All$Freq)
colnames(Race_Compare) <- c('Race','Women','50% of General population','Percentage Difference (%)')
Race_Compare

From the graph, we notices that majority of occupation were not meeting the criteria. To gain a better understand, we took into account he population size.

#Comparing with Histograms
DSOR_W <- data.frame(Data_sex_occp_race[which(Data_sex_occp_race$Gender==2),])
p <- ggplot(DSOR_W, aes(x=Race,y=MEAN, fill=factor(Race))) +
            geom_bar(position = "dodge",stat="identity")+
            facet_grid(. ~ Occupation)+
            geom_hline(aes(yintercept = MEAN),data=Data_sex_occp_women)+
            geom_text(aes(1,0.5,label = "Women's Mean"))+
            labs(title = "Histogram by Race of Females compared to Women's Mean")
ggplotly(p)

is.na() applied to non-(list or vector) of type 'NULL'is.na() applied to non-(list or vector) of type 'NULL'is.na() applied to non-(list or vector) of type 'NULL'is.na() applied to non-(list or vector) of type 'NULL'is.na() applied to non-(list or vector) of type 'NULL'is.na() applied to non-(list or vector) of type 'NULL'is.na() applied to non-(list or vector) of type 'NULL'is.na() applied to non-(list or vector) of type 'NULL'

Oberservations:

1. Asians women earn significantly higher income than other races.

2. The average for women is very close the White women because they account for most of the workforce.

3. There are no groups of women that are higher than the median of men

From the plot, we also see a lot of minority groups toward the bottom. The frequency for these froups were often very small. To help gain a better understanding, we took out all groups with less than 50 people in the category.

x <- ftable(Data_Women$Race,Data_Women$Occupation)
x <- as.data.frame(x)
x <- x[-31,] #frequency of zero
DSOR_W_Freq <- data.frame(DSOR_W,x$Freq)
Base <- ggplot(DSOR_W_Freq, aes(x=Race,y=MEAN),colour = factor(DSOR_W_Freq$Race)) +
            facet_grid(. ~ Occupation)+
            geom_point(aes(size = DSOR_W_Freq$x.Freq,color = factor(DSOR_W_Freq$Race)))+
            labs(title = "Histogram by Race of Females using Size")+
            geom_hline(aes(yintercept = MEAN),colour="deeppink1",data=Data_sex_occp_women)+
            geom_hline(aes(yintercept = MEAN),colour="blue",data=Data_sex_occp_men)+
            geom_hline(aes(yintercept = MEAN),colour="green",data=Data_occu)
ggplotly(Base)

Obeservations:

1. Except for Asians and Blacks, all other minority groups earned significantly less than the median.

2. Black Women earn less than White womens for all occupations, except legal.

3. Lawyers would have the most equal pay relative to the rest of the occupations.

From our analysis, we learned the following breakdown of women in the workforce:

Level	Race of Women
Below	Black, American Indian, Tribe specified, Native Hawaiian, Other, Two or more races
Meet	Asian, White
Exceed	None

Through these results, their is a significant gap in income for women regardless to race.

Appendix

Using Survey Weights

We were trying to determine whether or not to use survey weights? Would the data be significantly impacted?

DSOR_W_Freq_filtered <- filter(DSOR_W_Freq, x.Freq > 50)
Base_fill <- ggplot(DSOR_W_Freq_filtered, aes(x=Race,y=MEAN),colour = factor(DSOR_W_Freq_filtered$Race)) +
            facet_grid(. ~ Occupation)+
            geom_point(aes(size = DSOR_W_Freq_filtered$x.Freq,color = factor(DSOR_W_Freq_filtered$Race)))+
            labs(title = "Histogram by Race of Females using Size")+
            geom_hline(aes(yintercept = MEAN),colour="deeppink1",data=Data_sex_occp_women)+
            geom_hline(aes(yintercept = MEAN),colour="blue",data=Data_sex_occp_men)+
            geom_hline(aes(yintercept = MEAN),colour="green",data=Data_occu)
ggplotly(Base_fill)

From our results, the survey weights would be insignificant relative to the data. The difference for each income level between unweighted and weighted is too small.